Long-range correlations and coherent structures in magnetohydrodynamic equilibria 
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The equilibrium theory of the 2D magnetohydrodynamic equations is derived, accounting for the 
full infinite hierarchies of conserved integrals. An exact description in terms of two coupled elastic 
membranes emerges, producing long-ranged correlations between the magnetic and velocity fields. 
This is quite different from the results of previous variational treatments, which relied on a local 
product ansatz for the thermodynamic Gibbs distribution. The equilibria display the same type of 
coherent structures, such as compact eddies and zonal jets, previously found in pure fluid equilibria. 
Possible consequences of this for recent simulations of the solar tachocline are discussed. 

PACS numbers: 47.10.-g, 05.70.Ln, 05.90.+m, 52.30.-q 



The ideal magnetohydrodynamic (MHD) equations 

<9 t v + (v • V)v = -Vp + J x B 

d t B = V x (v x B) (1) 

describe the evolution the velocity field v of a perfectly 
conducting fluid under the influence of pressure gradients 
and the magnetic Lorentz force, and advection of the 
magnetic field B by the velocity field The equations 
are closed through Ampere's law J = V x B and the 
incompressibility constraints V • v = 0, V • B = 0. 

A number of systems such as the solar tachocline Q , 
and possibly the Earth's core- mantle boundary are 
approximately governed by a 2D approximation in which 
v, B are horizontal, depending only on the horizontal co- 
ordinates r = (x, y) in a domain D of the icy-plane (with 
ID boundary dD) [J]. The current J = Jz and vorticity 
V x v = wz may then be treated as scalars. One may 
also express B,v in terms of a potential A and stream 
function -0: B = Vx (Az), v = V x (ipz), with the 2D 
Laplacian relationships J = — V 2 A, us = —V 2 ip. The 
equations then reduce to two scalar equations 



dtuj + v- V(w + /) = B-VJ 
d t A + wVA = 



(2) 



in which the Coriolis parameter /(r) = 2£lsin(< j c), with 
latitude tp(y), accounts frame of reference rotating with 
angular velocity fl. In the beta-plane approximation one 
linearizes f(y) — fo + fty about some reference ifQ. The 
total energy is given by 



E = 



(|v| 2 + |B| 2 ) = i J d\(\V^ + \VAf 



(3) 

Its conservation requires lossless boundary conditions, 
e.g., periodic or "free-slip": B • n = 0, v • n = 0, with n 
the local normal to dD. The latter imply constant values 
of A, t/j on each connected component of dD. 

The pure 2D fluid equation (B = 0) exhibits a tur- 
bulent inverse cascade: small scale vortices self-organize 
into large scale flows, limited only by the domain size. 



In a forced system, this leads to growing flows, limited 
only by dissipation. These take the form of strong jet- 
like structures, or large coherent vortices, which have 
strong implications for geophysical flow stability and 
global transport. The weather bands and Great Red 
Spot of Jupiter are famous examples. For weak driv- 
ing and dissipation, these flows may be modeled as near- 
equilibrium. Application of statistical mechanics to the 
fluid equations indeed produces such structures [5rfl2|. 

Similar flows in the solar tachocline, which sharply di- 
vides the rigidly rotating interior radiation zone from the 
differentially rotating outer convection zone, would have 
strong implications for angular momentum transport be- 
tween the two zones [2]. Recent 2D MHD simulations 
d, however, have found that the presence of even 
weak B tends to break up large scale flows. Although 
B indeed destroys the infinite set of vorticity conserva- 
tion laws that lead to the fluid inverse cascade, they are 
replaced by an entirely new infinite set of magnetic con- 
servation laws 13[. This motivates the present investi- 



gation of the equilibria allowed by these new laws, and 
what conditions might limit their amplitude. 

It will be shown that the new conservation laws in- 
deed produce large scale flows, but of distinctly different 
character. In contrast to Euler flow, where all initial en- 
ergy cascades to large scales, in 2D MHD a finite fraction 
flows to small scales, generating large microscale fluctua- 
tions. For small initial B, these can be much larger than 
the mean flow 0, 0]- However, dissipation effects differ- 
entially tend to suppress these, and some form of driv- 
ing might amplify the mean flow to more visible levels. 
The results here may then point to more careful stud- 
ies of the types driving and initial conditions that could 
produce stronger large scale flows, and their physical re- 
adability. In addition, unlike Euler flow, in which the 
small scale fluctuations are completely uncorrelated on 
the microscale [5, 7] , the 2D MHD flows are shown to ex- 
hibit long-ranged power-law correlations, that could also 
be confirmed in simulations. Such correlations were in- 
visible to earlier treatments which used a loc al p roduct 
approximation for the Gibbs distribution 15|, [l6[ . 
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The MHD equations are nonlinear, thereby produc- 
ing turbulent initial evolution from a spatially complex 
initial condition. However, at late time, the flow may 
equilibrate to a steady state. Unusually, such 2D fluid 
steady states are not necessarily quiescent, but can ex- 
hibit spatial structure such as large scale vortices or zonal 
jets (HHl^ . The origin, and diversity, of these structures 
lies in the infinite number of conservation laws, beyond 
the usual energy (along with momentum or angular mo- 
mentum if the relevant translation or rotation symmetry 
exists) constraining the flow [l3j]. If one sets B = 0, 
the first line of ([2]) is the Euler equation, and the re- 
sulting advective conservation of the potential vorticity 
Up = uj + / implies conservation of all spatial integrals 
(Casimirs) K g = J D d 2 g(ui p )d 2 r, with g an arbitrary ID 
function. However, any nonzero B • V J breaks this con- 
servation replacing it by advective conservation of 
A [second line of ((5}]. The new Casimirs are [l3| 



Jo 



d 2 rg(A), K g 



d 2 r(u + f)g(A) (4) 



D 



for arbitrary ID functions g. Integrating by parts, one 
may also replace ujg(A) — > (v • Ti)g'(A). These may be 
parameterized as conservation of the functions j(cr), k(a) 
for all a obtained using g(s) = 5(s — a). Finally, any 
conserved momentum may be expressed in the form 



P 



d 2 rauj = / d 2 r(Va) x v 



d 2 rS7a-Vip, (5) 



for some fixed function a(r) [17 



Under the usual ergodic assumptions (whose validity 
is far from obvious, and known to be violated for some 
initial conditions |18j|). the equilibrium statistics are ob- 
tained from the grand canonical partition function 



Z\J3,\,n,v] = / D[A,i)]e-W^\ 



(6) 

and associated free energy F = — T\n(Z). The functional 
integral is over all fields A,ip (obeying the appropriate 
boundary conditions), /3 — 1/T is the inverse tempera- 
ture, and the Gibbs functional is Q = E — J M — K v — XP: 

G = J d 2 rl^\VA\ 2 + ^\Vip\ 2 -i/(A)VA-Vip 

+ \Va>Vil>-\n(A) + MA)]\, (7) 



where Lagrange multipliers have been introduced for each 
conserved integral: functions /i(s), v(s) and parameters 
/3, A are adjusted to obtain the values of j(a), k(a), P, E 
defined by the initial flow. The physical model associated 
with g is that of two membranes, with "heights" A,ip, 
and unit surface tension, coupled through their gradients. 
The P term also acts to bias Vi/>; in particular, g favors 
B parallel to v'(A)\ + AV x a. The term fj,(A) + fv{A) 



is an external potential, confining A near its minimum, 
and depends (smoothly) on position through /(r). 

The equilibria governed by g have been previously in- 
vestigated 15|, LL6| using a variational approach in which 
v, B were treated as spatially uncorrelated. The exact 
physics of C/, however, dictates something quite differ- 
ent. The elastic interactions generate very long-ranged 
correlations, with distortion in the surfaces interacting in 
a Coulomb-like fashion, leading to log-divergent fluctua- 
tions, and dipole-like correlations for v,B. 

To understand the nature of the equilibrium states, it is 
extremely useful begin with a discrete spatial mesh, and 
then consider the continuum limit. To this end, using a 
square lattice with mesh size a one obtains 

PSa = ^£{^-^') 2 + (^-^)) 2 

~ y^ + v'iAMAi-A^i-^) 

+ 2A(a 4 - Oj)(t/>t -ipj)} 

- 0a 2 Yy(Ai) + fiv{Ai)], (8) 

i 

in which (i, j) are nearest neighbors. The functional 
integral D[A,ip] — > J\ i f dAidipi now becomes (up to 
overall normalization) an independent product [19(. We 



-\A 



i+y 



Aj , Aj 



define also B.; = 
a~ 1 (ipi+y ~ ipi,ipi — i>i+z)- It is apparent here that if 
P remains finite as a — > then the first term will yield 
0(Vt) fluctuations between neighboring sites, and the 
second term becomes negligible [20( . Well defined hydro- 
dynamic equilibria, with nontrivial competition between 
kinetic and potential energies requires that /3 scale with 
a. Specifically, taking /3(a) = f3/a 2 yields 0{aT 1 / 2 ) site- 
to-site fluctuations, and continuous A(r),ip(r). On the 
other hand B^, v, will have finite, 0(Vt), site-to-site 
fluctuations, yielding non-differentiable continuum A,ip. 
With this scaling there are two surviving contributions to 
g as a — > 0. Decomposing A = Aq + SA, ip = ?/>o + Sip in 
which Aq = (A), tpo = (tp) are the equilibrium averages 
(self-consistently determined below) one obtains 



Q[A,ip] = g[A ,ip ] + g 2 [6A,8t,A ] 
[|VM| 2 



0(a), 



g 2 = » d 



- 2 l /(A )(V<L4) • (V61P)] 
d 2 r[[l + v'{A Q )]\V5 ( p-\ 2 
+ [l-i,'(A )]|V^+| 2 ], 



(9) 



in which Sdr^ = (SA ± 8ip)/y/2 are independent Gaus- 
sian fields [2l|. We have reverted, for compactness, to 
continuum notation. For smooth Aq, ipQ, all other terms, 
including those linear in SA, Sip, vanish with a — ¥ 0. From 
C/2 one therefore obtains the exact free energy 



(10) 



3 



in which the Gaussian fluctuation free energy Ti has a 
well defined a — > limit. It is not computable in closed 
form for general A , but for constant v'{Aq) = v' Q 
periodic boundary conditions one obtains 



Current density J (x,y); y = -50, fi = 0.05 



\D\ 



BZ 



d 2 k 



• 111 



L 

2tt 



1 - vgEQt) 



E(k) = 4[sm 2 (k x /2) + sm 2 (k y /2)] 



(11) 



in which \D\ is the area of D, and the Brillouin zone is 
defined by tt < k x , k y < tt. More generally, Ti is obtained 
from the log-sum of the eigenvalues of the (generalized 
Laplacian) operators ± = — V • [1 =F v'{Aq)]V . 

The divergence of j3 — j3 /a 2 in (|6|) means that ^0,^0 
are determined self-consistently by minimizing J- . The 
conditions SQo / Stpo(r) = 0, 5(Q + J r 2 )/M (r) = pro- 
duce, respectively, the equilibrium equations 



v = u'(A )B - AV x a 
u + v'(A )J = ii'(Ao) + fi/(Ao) 

+ u"(A )j(r,r;A ), 



(12) 



where Bo = V x (Aoz), etc., are the other derived equilib- 
rium fields. The first equation provides a direct relation 
between the equilibrium velocity and magnetic field: Bo 
is colinear with vo up to a mean flow subtraction. Sub- 
stituting the curl of the first equation into the second, 
one obtains a closed equation for Aq: 



-V x {[1 - is'(A ) 2 ]B Q } - fi'(A ) + fu'{A ) + XV 2 a 
+ i[^(Ao) 2 ]'|Bo| 2 + ^(A )7(r,r). (13) 



The fluctuation-derived nonlocal term at the end is ob- 
tained from the cross-correlation function 

7 (r,r') = (W5A(r) ■ V^(r')) 2 - (<JB(r) • 6v(r')) 2 . (14) 

If one defines the Green functions of the operators O^, 



- V • [1 T ^'(A )]VG ± (r, r') = S(r - r'), 
which characterize the membrane fluctuations, 



G ± (r,r') 



r ([^ ± (r)-^ ± (r')] 2 ) 2 , 



then one obtains the relation 



1. 



7 (r,r') = -V-V'[G+(r,r')-G-(r,r')]. 



(15) 



(16) 



(17) 



Solutions to ([T5|) yield the electrostatic potential gener- 
ated by a unit charge at r' with spatially varying di- 
electric function e ± (r) = 1 =p v'(A ). By Gauss's law, 
the displacement field e ± VG ± has unit flux through any 
bounding contour, and it follows that G ± (r — r') ~ 
(l/2ne ± ) ln(|r — r'\/a) grows logarithmically with dis- 
tance. This unbounded wandering reflects the usual ther- 
mal roughening result for 2D membranes [22J. 
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FIG. 1: Left: Equilibrium jet solution obtained using peri- 
odic boundary conditions in x. Right: Equilibrium vortex 
solution obtained using free slip boundary conditions in both 
x and y. Vortex solutions exist as well for periodic bound- 
ary conditions, but the numerics are more involved. Stream 
lines and magnetic field lines are level curves of constant po- 
tential Aq (bottom panels) or current density Jo (top pan- 
els). For simplicity we set T = (no microscale fluctua- 
tions), / = (no rotation), A = (no net mean flow), take 
constant v'(A) = p' (hence Bo = vo/fo), and use the form 
H{A) = - 7 _1 ln[l + e- 7(A_All) ]. For large 7 < 0, the latter 
tends to produce a two level system, with Jo = ujo/ v' q taking 
values or q = (1 — f 2 ) -1 , and the parameter fj,i controls 
the proportion of each level. Plotted are actually the scaled 
quantities Ao/q, Jo/q, with parameters 70 = 79, no = fiiq 
specified in the titles. The results are then independent of q 
(and hence of the choice of Uq). 



For a translation invariant system (i.e., periodic, uni- 
form), one may replace V' — > —V and (|15|) implies that 
7 = for r ^ r', and 7 (r,r) = v' /2j3[l - vg]. For this 
case, the locality assumption [l5[ is valid [23}. Otherwise 
one obtains nonlocal contributions to 7(1*, r), even for a 
uniform system with free slip boundary conditions: non- 
local contributions then arise from image charges outside 
D that enforce the boundary conditions on G^. 

Once the equilibrium solution is obtained, the values of 
the conserved variables may be obtained from derivatives 
of T at fixed Aq, ^>o : 



E 



1 



D 



d 2 r[|v | 2 + |B | 2 + £ (r,r)] 



P = -—=J D drau Jo 



3W 



ST 
ST 



d 2 rS[a - A (r)} 



k ^ = -Jrh= / d 2 r{(u + f)5[<j~A Q (T)} 

- 7 (r,r)5'[a-A (r)]}. (18) 

where e(r, r') = — V • V'[G + (r, r') + G~ (r, r')] determines 
the microscale fluctuation part of the energy. The ex- 
ternum condition ensures that all derivatives with re- 
spect to the implicit dependence of Ao,ij)Q on the La- 
grange multipliers cancels out. 
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Comparing Q , it is seen, due to continuity of A, that 
j'(cr) is a large-scale quantity: it is its own equilibrium 
average. Hence the level sets of the initial A(r,t = 0)) 
are exactly preserved in the equilibrium A (r). One im- 
plication is that a low-amplitude initial A will produce an 
identically low-amplitude Aq. However, e.g., for a spa- 
tially irregular initial condition with comparatively large 
energy, one will have e 3> |v | 2 ,B | 2 and the physical 
equilibrium fields will be masked by fluctuations. This 
may help explain what is observed in simulations 0, H| ■ 
For increasing initial A, especially in the presence of forc- 
ing and weak dissipation, one may expect clear equilib- 
ria to emerge more strongly, but the exact conditions 
required for this remain to be determined. 

Equations (fl2|) . (fT7| and (fTSj) are the basic results of 
this paper. These equations are nonlocal and highly non- 
linear, and generally require a numerical solution. How- 
ever, some general properties may be inferred from the 
basic form of the effective Hamiltonian (|9|). The Qq term 
reflects a straightforward classical surface tension min- 
imization problem. The fluctuations locally stretch the 
membranes and hence act to renormalize the surface ten- 
sion, but in a way the depends self-consistently on the av- 
erage membrane position Ao,ipo, which in turn respond 
to the external forces provided by the fi+fv and Xa terms 
in ([?]) ■ Thus, a stretched region of a membrane may be 
expected to have reduced amplitude fluctuations. The 
minimization of T accounts fully for both effects. Figure 
[1] shows examples of very simple equilibrium solutions 
in a rectangular domain (with model parameters speci- 
fied in the caption) , demonstrating the existence of both 
vortex and jet structures. 

Unlike the Euler case, where the conserved integrals 
ensure bounded to = V x v, thereby generating continu- 
ous v and differentiable i/>, for 2D MHD only continuity 
of A, ip are provided, and this occurs now in response to 
the squared-gradients in the energy, not in response to 
the conserved integrals. The latter also ensures bounded 
<5v,<5B, but this in turn allows for finite microscale en- 
ergy density e. The second derivatives uj, J then have 
unbounded fluctuations. Since simulations often propa- 
gate A,lu using ([2]), extra care may then be required to 
ensure reasonable equilibration. In particular, if (e.g., 
hyperviscous) dissipation acts too strongly to quell the 
micro-fluctuations, it may also bleed energy out of the 
large scale flow. There could also be physical analogues 
of this effect, depending on the precise nature, e.g., of 
the true solar tachocline dissipation mechanisms. These, 
and presumably many other, considerations must enter 
the implications of the theory developed here. 
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